Impurity probe of topological superfluid in one-dimensional spin-orbit coupled atomic 
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We investigate theoretically non-magnetic impurity scattering in a one-dimensional atomic topo- 
logical superfluid in harmonic traps, by solving self-consistently the microscopic Bogoliubov-de 
Gennes equation. In sharp contrast to topologically trivial Bardeen-Cooper-Schrieffer s-wave su- 
perfluid, topological superfluid can host a mid-gap state that is bound to localized non-magnetic 
impurity. For strong impurity scattering, the bound state becomes universal, with nearly zero energy 
and a wave-function that closely follows the symmetry of that of Majorana fermions. We propose 
that the observation of such a universal bound state could be a useful evidence for characterizing the 
topolgoical nature of topological superfluids. Our prediction is applicable to an ultracold resonantly- 
interacting Fermi gas of ' 
lattice. 



K atoms with spin-orbit coupling confined in a two-dimensional optical 
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I. INTRODUCTION 

Impurity scattering plays an important role in under- 
standing the quantum state of hosting systems [l| . This 
is particularly significant in solid state systems, where 
impurity scattering and disorder are intrinsic. In super- 
conductors, the study of impurity effects has the poten- 
tial to uncover the nature and origin of the supercon- 
ducting state [2]. In strongly correlated electronic sys- 
tems near quantum critical points, where several types 
of ordering compete in a delicate balance, the study of 
impurity scatterings has the power to underpin in favor 
of one of the orders [3]. In this work, we aim to investi- 
gate theoretically impurity scattering in one-dimensional 
(ID) topological superfluids. We show that an impurity- 
induced bound state will provide a sensitive probe for the 
topological order in such systems. 

Topological superfluid is a novel state of quantum mat- 
ter [J| , which is gapped in the bulk, but hosts non-trivial 
zero-energy surface states - the called Majorana fermions 
IH - near its boundary. It has attracted great at- 
tentions in recent years because of its potential appli- 
cation in topological quantum computation and quan- 
tum information [7j, [8j. The realization of topological 
superfluids and the manipulation of Majorana fermions 
are currently the most hot research topic in a variety 
fields of physics, ranging from condensed matter physics 
to ultracold atomic systems. Till now, indirect evi- 
dence of the existence of topological superfluids in hybrid 
superconductor-semiconductor InSb or InAs nanowires 
has been reported [9l-[TT| . Theoretical schemes of process- 
ing topological quantum information in such nanowire 
devices have also been proposed (l2l. Il3j. 
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Our investigation of impurity scattering in ID topo- 
logical superfluids is strongly motivated by the rapid ex- 
perimental progress On one hand, impurity scat- 
tering is un-avoidable in InSb or InAs nanowires. A real- 
istic simulation of impurity scattering may therefore be 
useful for future solid-state experiments. On the other 
hand, we anticipate that impurity may induce new ex- 
otic bound state, thus providing a clear local probe of 
the topological nature of the systems that we consider. 

In this paper, we use a ID spin-orbit coupled atomic 
Fermi gas to model ID topological superfluids (l^ - tl6j . in- 
stead of considering nanowire devices used in solid-state 
[gl-TlTjl . This is because we have unprecedented control- 
lability with ultracold atomic gases [l7j . By using mag- 
netic Feshbach resonances, the interatomic interactions 
can be precisely tuned (la ]. Using the technique of op- 
tical lattices, artificial ID and 2D environments can be 
easily created pj| |2(j. The spin-orbit coupling, which 
is the necessary ingredient of a realistic topological su- 
perfluid, can also be engineered with arbitrary strength 
[2lll22|. Thus, ultracold spin-orbit coupled atomic Fermi 
gas is arguably the best candidate to simulate the desired 
topological superfluids. Furthermore, even though cold 
atom systems are intrinsically clean, individual impuri- 
ties can be realized using off-resonant dimple laser light 
or another species of atoms or ions [23j • The disorder ef- 
fect of many randomly distributed impurities can also be 
created by employing quasiperiodic bichromatic lattices 
or laser speckles |24| . 

We investigate the impurity effect in ID spin-orbit 
coupled atomic Fermi gas of 40 K atoms by solving 
self-consistently the microscopic Bogoliubov-de Gennes 
(BdG) equation, with realistic experimental parameters. 
We observe the existence of mid-gap state that is bound 
to localized non-magnetic impurity. For strong impurity 
scattering, the bound state tends to be universal, with 
nearly zero energy and a wave-function that closely fol- 
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lows the symmetry of that of Majorana fermions. This 
feature is clearly absent in topologically trivial superfiu- 
ids. Therefore, we argue that the observation of such 
a universal bound state would be a useful evidence for 
characterizing the topological nature of topological su- 
perfluids. We note that, mid-gap bound state induced 
by non-magnetic impurity has also been predicted in ID 
spin-orbit coupled superconductors, by using non-self- 
consistent T-matrix theory [25]. The effect of magnetic 
impurity in 2D spin-orbit coupled Fermi gases has also 
been studied analytically using T-matrix formalism [2(| . 

Our paper is arranged as follows. In the next section 
(Sec. II), we introduce briefly the model Hamiltonian 
and the solution of BdG equations, and then present a 
phase diagram for a given set of experimental parame- 
ters. In Sec. Ill, we study non-magnetic impurity scat- 
terings and show the emergence of universal bound state 
in the strong scattering limit. The properties of such a 
universal bound state are analyzed in greater detail. To 
better simulate the realistic experimental setup, we also 
consider an extended impurity with gaussian-shape scat- 
tering potential. Finally, we summarize in Sec. IV. The 
detailed numerical procedure of solving BdG equations is 
listed in the Appendix A, together with a careful check 
on numerical accuracy. 



II. MODEL HAMILTONIAN AND BDG 
EQUATIONS 

The framework of our theoretical approach has been 
briefly described in our previous work [15| . Here, we em- 
phasize on the experimental origin of the model Hamil- 
tonian and generalize the theoretical approach to include 
a classical non-magnetic impurity. A detailed discussion 
on the numerical procedure is given in the Appendix A. 



A. ID spin-orbit coupled Fermi gas 

Let us consider a spin-orbit-coupled Fermi gas of 40 K 
atoms in harmonic traps, realized recently at Shanxi Uni- 
versity [2l|. We assume additional confinement due to a 
very deep 2D optical lattice in the transverse y — z plane, 
which restricts the motion of atoms to the a;-axis. The 
spin-orbit coupling is created by two counter propagat- 
ing Raman laser beams that couple the two spin states 
of the system along the x-axis [2l[ . Near the Feshbach 
resonance Bo ~ 202.20 G, the quasi-lD Fermi system 
may be described by a single-channel model Hamiltonian 
H = H Q + H int , where 
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(1) 



is the single-particle Hamiltonian in the presence of Ra- 
man process and 



Hi, 



gin 



dxM>\ (x) *{ (x) j. (x) * t (x) (2) 



is the interaction Hamiltonian describing the contact in- 
teraction between two spin states. Here, the pseudospins 
(7 =t, I denote the two hyperfme states, and fy a (x) is 
the Fermi field operator that annihilates an atom with 
mass m at position x in the spin a state. The chemical 
potential /i is determined by the total number of atoms 
N in the system. For the two-photon Raman process, flu 
is the coupling strength of Raman beams, kji — 2tt / \r 
is determined by the wave length of two lasers and 
therefore 2hkji is the momentum transfer during the pro- 
cess. The trapping potential Vr (x) = mu> 2 x 2 /2 refers to 
the harmonic trap with an oscillation frequency lj — uj x 
in the axial direction. In such a quasi-one dimensional 
geometry, it is shown by Bergeman et al. that the 
scattering properties of the atoms can be well described 
using a contact potential guj5(x), where the ID effective 
coupling constant g±o < may be expressed through the 
3D scattering length 0315, 



9id = 



2h 2 a 3D 



1 



(1 - Aa 3D /aj_) ' 



(3) 



where a± = ^h/(muj±) is the characteristic oscillator 
length in the transverse axis, for a given transverse trap- 
ping frequency ujj_ set by the deep 2D optical lattice. 
The constant A = -C(1/2)/a/2 ~ 1.0326 is respon- 
sible for the confinement induced Feshbach resonance 
[27l | . which changes the scattering properties dramati- 
cally when the 3D scattering length is comparable to the 
transverse oscillator length. It is also convenient to ex- 
press gin in terms of an effectivelD scattering length, 
giD = -2ft 2 / (mai D ), where a 1D = -(a 2 L /a 3D )(l - 
Aa,3E>/a±) > 0. The interatomic interaction can then 
be described by a dimensionless interaction parameter 
7 = a/[-K\Na\D\, where a = \J h/{muj) is the oscillator 
length in the x-axis. Near the Feshbach resonance, the 
typical value of the interaction parameter 7 is about 5 

[IE BUI. 

To illustrate how the spin-orbit coupling is induced by 
the two-photon Raman process, it is useful to remove 
the spatial dependence of the Raman coupling term, by 
taking the following local gauge transformation, 



* t (as) 
*4 (*) 



_ e +ik R x 



lkRX ^i (x) . 



(4) 
(5) 



Using the new field operators ifo (x) and ip± (x), we can 
recast the single-particle Hamiltonian as 



Ho = 



Ho 



dx 



$f ( x ) ' % ( x ) 



He 
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where we have absorbed a constant energy shift En = 
f?k 2 R /(2m) (the recoil energy) in the chemical poten- 
tial n, and have defined the momentum operator k x = 
—id/dx, the spin-orbit coupling constant A = H 2 kn/m 
and an effective Zeeman field h = flji/2. a x and a z are 
Pauli's matrices. The spin-orbit coupling in the Hamil- 
tonian can be regarded as an equal-weight combina- 
tion of Rashba and Dresselhaus spin-orbit coupling (i.e., 
Xk x a y ). This is evident after we take the second local 
gauge transformation, 



ipf (x) 



1 

71 
i 

71 



bPt ( x ) - % ( x )} . 

[tpf (x) + iip± (x)] , 



(8) 
(9) 



with which the single-particle Hamiltonian becomes, 



Hn = dx 



rtn 



(10) 



Ho 



h 2 d 2 

2m dx 2 



+ Vt (x) — fi — ha z + Xk x a y . (11) 



The form of the interaction Hamiltonian is invariant after 
two gauge transformations, i.e., 

H in t = gm / dxip\ (x) ip\ (x) ipi (x) (x) . (12) 



We note that the operator of total density n(x) 
J2 a ( x ) ( x ) = J2a Va ( x ) V'cr {x) is also invariant in 
the gauge transformation. 



B. Impurity scattering Hamiltonian 

Now we add the non-magnetic impurity scattering 
term, 

H imp = / dxVimp (x) Ipl (x) lj} a (x) , (13) 

to the total Hamiltonian. The non-magnetic scattering 
can be realized experimentally by using an off-resonant 
dimple laser light. We consider either a localized scatter- 
ing potential at position Xq, 



V imp (x) = V imp S (X - Xq) , 



(14) 



or an extend potential with a width d in the gaussian 
line-shape, 



Vimp (^0 



2nd 



exp[ 



(x - x y 
2d 2 



(15) 



The strength of the impurity scattering is given by Vi mp . 
In the narrow width limit d — > 0, the gaussian potential 
returns back to the delta-like potential. We may place 
the impurity at arbitrary position, as long as the Fermi 



system is locally in the topological superfluid state. To 
be concrete, we shall set xo = 0. 

We may also consider a magnetic impurity scatter- 
ing in the form, H lmp = J dxVi mp (x) [V>j (x) ipf (x) - 

ip\ (x) (x)]. However, it is of theoretical interest only. 
The field operator of density difference is not invariant 
in the second local gauge transformation. Thus, exper- 
imentally the magnetic impurity scattering potential is 
more difficult to realize. 



C. Bogoliubov-de Gennes equation 

We use the standard mean-field theory to solve the 
model Hamiltonian. By introducing a real order param- 
eter A (a;) = —giD (V'J, ( x ) the interaction Hamil- 
tonian is decoupled as, 



Hint 



dx 



A (x) ipfyl (x) + H.c. + 



|A(x)r 

9lD 



(16) 

It is then convenient to introduce a Nambu spinor i/>(x) = 
[V>t (x) , tpi (x) , iftt (x) , ip] (x)] T and rewrite the mean- 
field Hamiltonian in a compact form, 



H mf dx^U BdG ip(x) + TtHs - J dx 



\A(x) 



9lD 



(17) 



_ where 



H 
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H s -h -Xd/dx 
Xd/dx H s + h A(x) 

A*(x) -Hs + h 

-A*(x) -Xd/dx 



-A(x) 


Xd/dx 

-u s -h 



and 



(18) 



H s (x) 



h 2 d 2 
2m dx 2 



+ — wV + V lmp {x). (19) 



The term TrHs m H rn f results from the anti- 
commutativity of Fermi field operators. 

The mean-field Hamiltonian Eq. (TTT)) can be diago- 
nalized by the standard Bogoliubov transformation. By 
defining the field operators a v for Bogoliubov quasipar- 
ticles, 

dx ^ [ U >Jri (x) VV (x) + v ar] (x) -4>l (x)] , (20) 
a 

we obtain that, 

|A(x)| 2 



I r 
H m f = - 2J E n u\u n + TtHs - / dx- 



9id 



(21) 



Here, $ v (x) = [u^ (x) , u iv (x) , u t „ (x) , v iv (x)] T and Ej, 
are respectively the wave-function and energy of Bogoli- 
ubov quasiparticles, satisfying the BdG equation, 



UsdG^n (x) = E n % (x) 



(22) 
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The BdG Hamiltonian Eq. (fT5|) includes the pair- 
ing gap function A (x) that should be determined self- 
consistently. For this purpose, we take the inverse Bo- 
goliubov transformation and obtain 

Vv (») = [ u *v ( x ) a v + v *v ( x ) a v\ ■ ( 23 ) 
n 

The gap function A (x) is then given by, 

(x) 4, (i) /(-£„)] , (24) 

where / (E) = l/[e E / kBT + 1] is the Fermi distribution 
function at temperature T. Accordingly, the total density 
take the form, 

n (*) = \ E [i^ ^)i 2 / + K (*)i 2 / (~ E v) ■ 

The chemical potential fi can be determined using the 
number equation, N — J dxn (x) . 

It is important to note that, the use of Nambu 
spinor representation enlarges the Hilbert space of the 
system. As a result, there is an intrinsic particle- 
hole symmetry in the Bogoliubov solutions [2^, |29| : 
for any "particle" solution with the wave-function 

$^ p) (x) = [tpf n (x) , ui n (x) , Vfr, (x) , v iv (x)] T and en- 
(p) 

orgy we can always find another part- 

ner "hole" solution with the wave-function <&^\x) = 

\ v \ v ( x ) > v t v ( x ) > u f v ( x ) > u ln ( X )] T and ener gy E v = 
—E^ < 0. In general, these two solutions correspond 
to the same physical state. To remove this redundancy, 
we have added an extra factor of 1/2 in the expressions 
for pairing gap function Eq. (|24[) and total density Eq. 
(HO). 

The Bogoliubov equation Eq. (fT5|) can be solved itera- 
tively with Eqs. (|2"4")l and by using a basis expansion 
method, together with a hybrid strategy that takes care 
of the high- lying energy states [lij |28|, [29J . A detailed dis- 
cussion on the numerical procedure and a self-consistent 
check on the numerical accuracy are outlined in the Ap- 
pendix A. 

D. Phase diagram in the absence of impurity 

In our previous study [15| , we have discussed the phase 
diagram of a weakly interacting spin-orbit coupled Fermi 
gas, with an interaction parameter 7 = tt/2 ~ 1.6. The 
real experiment, however, would be carried out near Fes- 
hbach resonances, where the typical interaction param- 
eter is 7 = 3 ~ 5 [20L l28l [29| . In this work, we take a 
realistic interaction parameter 7 = it ~ 3.2, despite the 
fact that our mean-field treatment would become less ac- 
curate. We consider a Fermi gas of N = 100 atoms in 




Figure 1: (color online) Phase diagram at T — (solid line) 
and T = 0.37> (gray dashed line), determined from the be- 
havior of the lowest energy in quasiparticle spectrum. With 
increasing the effective Zeeman field, the Fermi cloud changes 
from a standard BCS superfluid to a topologically non-trivial 
superfluid. The phase transition point is slightly affected by 
finite temperature. The inset shows the energy spectrum at 
h — 1.2Ef as a function of the position of quasiparticles. A 
zero-energy quasiparticle (i.e., Majorana fermion) at the trap 
edge has been highlighted by a big dark circle. Here, we char- 
acterize approximately the position of a quasiparticle by using 
its wave-function: (x 2 ) = J dxx 2 ^^[u 2 , (x) + v 2 (a;)]. 

a single tube formed by a tight 2D optical lattice, and 
take the Thomas-Fermi energy Ep = k^Tp = (N /2)hw 
and Thomas- Fermi radius xp = yNa as the units for 
energy and length, respectively. For the spin-orbit cou- 
pling, we use a dimensionless parameter Xkp/Ep = 1, 
where hp — \/2mEp is the Thomas-Fermi wavevector. 

Fig. 1 presents the phase diagram at these parame- 
ters and at two temperatures T — and T — 0.3Tp, 
showing the well-known topological phase transition at 
a critical effective Zeeman field h c ~ l.Q5Ep. The dif- 
ferent phase is characterized by the lowest energy of Bo- 
goliubov quasiparticles, min{|i?^|}. At a small Zeeman 
field h < h c , the system is a standard Bardeen-Cooper- 
Schrieffer (BCS) superfluid, with a fully gapped quasi- 
particle energy spectrum (i.e., min{|_E r) |} > 0). Once 
h > h c , however, topologically non-trivial phase emerges. 
Though the quasiparticle energy spectrum is still gapped 
in the bulk, gapless excitations - Majorana fermions - ap- 
pear at the edges [33|, leading to an exponentially small 
lowest energy in the spectrum. This is fairly evident in 
the inset, where we plot the energy spectrum as a func- 
tion of the position of quasiparticles. 

To determine the critical Zeeman field h Cl we note that 
for a ho mog eneous spin-orbit coupled Fermi gas, it is 
given by [30ll3l| 

h c = vV + A 2 . (26) 

In harmonic traps as we consider here, the critical Zee- 
man field becomes position dependent. The local critical 
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Figure 2: (color online) (a) The pairing gap distribution func- 
tion A (a;) (dot-dashed line) and the criterion for a local topo- 
logical phase h > h c (x) (solid line) at h = 1.2Ef- The shaded 
cross-hatching highlights the topologically non-trivial area, 
(b) The linear contour plot of the local density of state at 
h — \.2Ef- At each trap edge, a series of edge states, in- 
cluding the zero-energy Majorana fermion mode, are clearly 
visible. 



Zeeman field, calculated using h c (x) — \J fi 2 (x) + A 2 (x), 
with the local chemical potential n(x) — [i — muj 2 x 2 /2 
and the local pairing gap A(x), increases monotonically 
towards the trap edge [34(. The Fermi cloud at position 
x will locally be in a topological state if the Zeeman field 
h > h c (x). For the parameters given in the above, this 
first happens at h ~ 1.05Ef, for which the local phase 
at the trap center (x — 0) starts to become topologi- 
cally non-trivial. In Fig. 2(a), we show the local pairing 
gap A(x) and the criterion for a local topological state, 
h > h c (x), at the Zeeman field h = 1.2Ep- At this field, 
the topological area is extended to the edge of the trap, as 
highlighted by a shaded cross-hatching. The appearance 
of Majorana fermion modes may be probed by measur- 
ing the local density of state through spatially resolved 
radio- frequency (rf) spectroscopy |15l. |23|. In Fig. 2(b), 



we present the local density of state at h — 1.2Ep, 
p(x,uj) = [\u an {x)\ 2 6 (u - E v ) 

(7Tj 

+ \v an {x)\ 2 8{lo + E n ) 



(27) 



At each of the two trap edges, we observe a series of edge 
states, whose dispersion relation is approximately given 
by [la ]. E n = y^AE, where n = 0, 1,2, • • ■ is a non- 
negative integer and AE is a characteristic energy scale 
set by the trapping frequency and recoil energy. The 
Majorana fermion modes with zero energy E n = are 
clearly visible. 



III. UNIVERSAL IMPURITY-INDUCED 
BOUND STATE 

We are now ready to investigate how Bogoliubov quasi- 
particles are affected by a non-magnetic impurity. Here- 
after, we focus on the topological state at ft. = 1.2i?p. 
For a topologically trivial state at h < 1.05Ep, we have 
checked numerically that quasiparticles are essentially 
not affected by the non-magnetic impurity scattering. 
This is in accord with the well-known Anderson's the- 
orem that potential scattering impurities are not pair- 
breakers in s-wave superconductors P, |32| . 



A. Impurity-induced mid-gap state 

In Fig. 3, we report the density profile and pair- 
ing gap distribution in the presence of a strong non- 
magnetic impurity with scattering potential strength, 



Vi, 



-0.30xfEf. Both of them are completely de- 



pleted at the impurity site x = 0. Accordingly, we 
observe the appearance of a new mid-gap state that is 
bound to the impurity, as shown in the inset for the 
spatial distribution of Bogliubov quasiparticles. This is 
clearly seen when we compare the quasiparticle spectrum 
without and with the non- magnetic impurity, i.e., the in- 
set in Fig. 1 and Fig. 3, respectively. Away from the 
impurity site, the distribution of Bogoliubov quasiparti- 
cles is also disturbed by the impurity. However, the series 
of edge states at the trap edge seems to be very robust 
against the impurity scattering. 

In Fig. 4, we show the local density of state p(x,u). 
The mid-gap bound state can be easily identified in spa- 
tially resolved rf spectroscopy, which is a cold-atom ana- 
log of scanning tunneling microscopy (STM). If such a 
bound state exists, one would observe a strong rf-signal 
at around origin and zero energy, which decays expo- 
nentially in space and energy. The maximum rf-signal, 
however, is located slightly away from the origin, as the 
total density is completely depleted right at the impurity 
site. 

The existence of a mid-gap state in the topological su- 
perfluid phase is certainly not consistent with Anderson's 
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Figure 3: (color online) Density profile and pairing gap 
distribution in the presence of a strong attractive non- 
magnetic impurity with scattering potential strength, Vi mp = 
— O.SOxfEf- The inset shows the spatial distribution of Bo- 
goliubov quasiparticle energy spectrum. The mid-gap bound 
state near the impurity site x — is highlighted by big blue 
circles. 

theorem [32| for potential scattering in s-wave supercon- 
ductors. However, it can be understood from the com- 
bined effect of the spin-orbit coupling and effective Zee- 
man field. Beyond the critical Zeeman field h c , the Fermi 
cloud is actually a p-wave-like superfluid (see, for exam- 
ple, the discussion in Sec. IIA of Ref. [16]). This is 
also the underlying reason why the cloud is in a topo- 
logical state. For superfluids with a non-zero angular 
momentum order parameter, non-magnetic impurity is a 
pair-breaker and would lead to a mid-gap bound state. 



B. Universal mid-gap state 

An impurity-induced bound state is not a unique fea- 
ture of topological superfluids, as it can also exist in su- 
perfluids with even-parity angular momentum order pa- 
rameter, such as d-wave and g-wave superfluids. Here, 
however, we argue that the existence of a deep, universal 
in-gap bound state in the limit of strong impurity scat- 
tering would be a robust feature of topological superflu- 
ids. Despite of the details of impurity scattering (i.e., 
non-magnetic or magnetic impurity, positive or attrac- 
tive scattering potential), we would observe exactly the 
same bound state, when the impurity scattering strength 
is strong enough. This argument is based on the consid- 
eration that a strong impurity will always deplete the 
atoms at the impurity site and hence create a vacuum 
area that is topologically trivial. Thus, at the interface 
between the topologically non-trivial and trivial areas, 
we would observe a pair of Majorana edge states [25| - 
the precursor of the universal bound state. Ideally, the 




Figure 4: (color online) (a) Density of state for a topological 
superfluid (h = 1.2Ef), at x — 0, 0.05:ef, and O.Ixf (from 
bottom to top). For better illustration, the curves have been 
off-set. The magnitude of the local density of state at x = 
and O.Ixf has been enlarged by a factor of 10. (b) Linear 
contour plot of local density of state. The impurity induced 
bound state is clearly visible near x = and uj = 0. 



energy of the universal bound state will be zero. 

In Fig. 5, we plot the energy of the mid-gap bound 
state as a function of the impurity scattering strength at 
h = 1.2Ep- Indeed, when the absolute value of the scat- 
tering strength Vi mp is sufficiently large, the energy of the 
bound state converges to a single value, E ~ 0.113Ao, 
where Ao ~ 0.464Ei? is the pairing gap at the trap cen- 
ter in the absence of impurity (see Fig. 3). We have 
also checked the case with a magnetic impurity and have 
found the same bound state energy (not shown in the fig- 
ure). The same bound state energy, found under differ- 
ent type of strong impurities, is a clear indication of the 
emergence of a universal impurity-induced bound state. 
It would also be a unique feature of the existence of a 
topological superfluid. 

We note, however, that the bound state energy is not 
precisely zero as we may anticipate from the Majorana 
edge-state picture as mentioned in the above. This is due 
to the fact that a pair of zero-energy Majorana fermions, 
localized at the same position (i.e., impurity site), could 
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Figure 5: (color online) The dependence of the bound state 
energy on the impurity strength for a topological superfluid 
at h — \.2Ef- The solid and empty circles show the results 
for attractive and repulsive potential scattering, respectively. 
The dashed lines gives the bound state energy at the infinitely 
large impurity strength, E ~ ±0.113Ao, obtained by an ex- 
trapolation. Here, Ao — 0.464_Ef is the pairing gap at the 
trap center without impurity. 



interfere with each other, leading to a small energy split- 
ting whose magnitude would depend on the detailed con- 
figuration of the Fermi cloud. In Fig. 6, we present the 
wave-function of the universal impurity-induced bound 
state. Indeed, the wave-function of the universal bound 
state can be viewed as the bond and anti-bond super- 
position of the wave- functions of two Majorana fermions, 
which satisfy the symmetry of u a (x) = v* (x) or u a (x) = 
—v* (x), respectively. 

We note also that the mid-gap state induced by 
non-magnetic impurities in topological superconduct- 
ing nanowires has recently been predicted by Sau and 
Demler, based on a non-self-consistent T-matrix and 
Green function method [25]. By increasing the impurity 
strength, it was reported that the bound state energy sat- 
urates to zero-energy, instead of converging to a nonzero 
value. In addition, a shallow bound state was predicted 
in the non-topological superconducting phase with spin- 
orbit coupling. These predictions are different from our 
numerical results. We ascribe these discrepancies to the 
lack of self-consistency in the T-matrix approach. 



C. Realistic gaussian-shape impurity 
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Figure 6: (color online) Wave-function of the universal bound 
state with energy E ~ ±0.113Ao, for a topological superfluid 
at h — \.2Ef- The wave-function at ±E may be regarded as 
the bond and anti-bond superposition of two Majorana wave- 
functions, which satisfy respectively the symmetry u a (x) = 
v% (x) (on the left side with x < 0) and u a (x) — —v% (x) (on 
the right side with x > 0). 



We consider so far a delta-like impurity scattering po- 
tential. In real experiments, the non-magnetic impurity 
would be simulated by an off- resonant dimple laser light, 
which has a finite width in space. Thus, it is more reason- 
able to simulate the impurity by using a gaussian-shape 
scattering potential. 




Figure 7: (color online) Linear contour plot of density of state 
at h — 1.2Ef, for an attractive or a repulsive gaussian-shape 
impurity scattering potential. Here, we take d = O.Isf and 
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Fig. 7 reports the linear contour plot of local density 
of state at h = \.2Ef for a strong attractive (a) and 
repulsive (b) gaussian-shape impurity potential. With 
a finite width d = O.Ixf, we observe a series of bound 
states in the vicinity of the impurity site. The lowest- 
energy bound state is close to the universal bound state 
that we find earlier with a delta-like impurity potential. 

To give some realistic parameters, let us consider a 
spin-orbit coupled Fermi gas of 40 K atoms confined to a 
tight 2D optical lattice, with an axial trapping frequency 
uj = 2ir x 116 Hz [21(. By assuming the number of atoms 



N = 100 in each tube, the Fermi energy or temperature 
is about 300 nK. We may take kp ~ 2fc^ and a Ra- 
man strength ~ 10_Er, where En is the recoil energy. 
We may anticipate a topological superfluid at tempera- 
ture T < 10 nK. The typical size of the Fermi cloud is 
about xf — 15 /im. Thus, we may use an off-resonant 
dimple laser with width d ~ 1.5 \xm to simulate the non- 
magnetic impurity. The strength of the impurity can be 
easily tuned by controlling the strength of the dimple 
laser light. With these parameters, we may be able to 
observe the universal impurity-induced bound state dis- 
cussed in the above. 



IV. CONCLUSIONS 

In summary, we have argued that a strong non- 
magnetic impurity will induce a universal bound state 
in topological superfluids. This provides a unique fea- 
ture to characterize the long-sought topological superflu- 
ids. We have proposed a realistic setup to observe such a 
universal impurity-induced bound state in atomic topo- 
logical superfluids, which are to be realized in spin-orbit 
coupled Fermi gases of 40 K atoms. The necessary con- 
ditions, including the realization of spin-orbit coupling 
by two-photon Raman process, the achievement of one- 
dimensional confinement by optical lattice, and the simu- 
lation of non-magnetic impurities using off-resonant dim- 
ple laser light, are all within the current experimental 
reach. Therefore, we anticipate our proposal will be real- 
ized soon at Shanxi University in China (2l| or elsewhere. 
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Appendix A: Solving the BdG equation in one 
dimension 

We solve the BdG equation Eq. (|22|) by expand- 
ing the Bogoliugbov wavefunctions u a (x) and v a (x) 
in the basis of ID harmonic oscillators fa(x) = 
(l/^TT 1 / 2 2ij\)H J (x)e- x ^ 2 , 

M-l 

u* (*) = U 'i^i ( x ) ' ( A1 ) 

3=0 
M-l 

v a (a?) = J2 V °i<t>i ( x ) ' ( A2 ) 

3=0 

Here, Hj (x) is the j-th Hermite polynomial and, for con- 
venience, we have used the natural unit in harmonic 
traps, m — h = lu = 1, so that the oscillator length 
a = \Jhj{miS) = 1 and the oscillator energy huj = 1 . On 
such a basis, the BdG Hamiltonian Eq. (|18|) is converted 
to a AM x AM secular matrix, 



h8j 



-R 



H 



BdG 



Rij 



-A, 



hSi 





A,;, 



A, 



-H s 







(A3) 



-Ri 



-H s 



hSj 



where the matrix elements, 

(i + 1/2 - n) Sij 



Ra 



V 13 

imp f 



(A4) 
(A5) 



To calculate efficiently the matrix elements 



f*™ dxfa(x)V imp (x)fa(x) and A, 



imp 

J_ °° dxfa(x)A(x)fa (x), we discretize space (-L/2, L/2) 
into N gr id equally spaced points, where the simula- 
tion length L and the number of grid N gr id should 
be sufficiently large so that the basis function fa(x) 
(j = 0,..,M — 1) can be accurately sampled. At the 
number of atoms N = 100, typically we take M = 500, 
Ng r id = 6400, and L = 7Qy/h/(rruS). The gaussian im- 
purity potential Vi mp (x) and pairing gap function A (x) , 



as well as the total density n(x), will be stored as an 
array of length N gr id- We note that, for a delta- like im- 
purity we immediately have V-^ p = V tmp 4>i (x ) fa (x Q ). 
By diagonalizing the AM x AM secular matrix Eq. 
(1A3[) . we obtain the quasiparticle energy E n and the 
eigenvector and V^j (j = 0,..,M — 1). The latter 
gives the quasiparticle wave- function u ar] [x) and v art {x) . 
Note that, the eigenvector U%j and have to satisfy 

the condition Y.ajWZj? + ( v * 3 ) 2 ] = h due to the 
normalization of the quasiparticle wavefunctions, i.e., 



1. 



In the practical calculation, due to computational lim- 
itation, we have to use a finite expansion basis. This is 
controlled by the cut-off M for the number of ID har- 
monic oscillators. Furthermore, we must impose a high 
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energy cut-off E c for the quasiparticle energy levels. To 
make our result cut-off independent, we adopt a hybrid 
approach, in which we solve the discrete BdG equation 
for the energy levels below the high energy cut-off E c . 
While above E c , we use a semiclassical plane- wave ap- 
proximation for the wavefunctions that should work very 
well for high-lying energy levels. For simplicity, to take 
the semiclassical approximation we may neglect the spin- 



orbit coupling term Rij in the BdG Hamiltonian Eq. 
(|A3[) . In the end, for the pairing gap function and the 
total density, we shall use the semiclassical expressions 
listed in the Sec. IVC of Ref. [28[. To summarize briefly, 
the contributions of discrete low-lying energy levels (la- 
beled by an index "77") and continuous high- lying energy 
levels to the total density are given by, 



rid 



(A6) 



\E V \<E C a 



and 



n c (x) 

•mil, \ 

Hc+h E a -h , 

respectively. For the pairing gap function, we have 




+00 +00 s 



de 



£ /Ve 2 - A 2 (x) - 1 



(A7) 



A(x) = 



e// ( \ 
9lD W 



E E K») ^ V U W f ( E n) + u ^ ( x ) v lv ( x ) / (- E n)] » 



(A8) 
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where the effective interaction strength g^J (ar) is deter- 
mined by, 



1. Check on the numerical accuracy 



gfj (x) 9iD 4tt?j J ' '^e 2 -A 2 (x) 

E c — h 



de- 



(A9) 



The numerical procedure of solving the BdG equation 
is therefore as follows. For a given set of parameters (N, 
9id, h, A, and T), we (1) start with an initial guess or 
a previously determined better estimate for A(x), (2) 
solve Eq. (|A9I) for the effective coupling constant, (3) 
then solve Eq. (|A3[) for all the quasiparticle wavefunc- 
tions up to the chosen energy cut-off to find u arj (x) and 
v ar) (x) , and finally determine an improved value for the 
order parameter from Eq. (|A8[) . During the iteration, 
the total density n(x) — rid(%) + n c (x) is updated. The 
chemical potentials 11 is adjusted slightly in each itera- 
tive step to enforce the number-conservation condition 
f-00 dxn(x)—N, until final convergence is reached. 



We have checked carefully the numerical accuracy of 
our hybrid approach at different sets of parameters and 
at both zero temperature and finite temperatures. In 
Fig. 8, we check the dependence on the cut-off energy E c 
at h — 1.2Ep in the absence of impurity scattering. The 
pairing gap function becomes essentially independent on 
E c once E c > 6Ep, with a relative error less than 1%. 
The cut-off energy dependence for the total density is 
even weaker (not shown in the figure). Thus, we con- 
clude that our hybrid calculation is quantitatively reliable 
with E c = 6Ep. At this energy cut-off, each iteration in 
the self-consistent calculation takes approximately sev- 
eral minutes, by using a standard desktop computer. The 
convergence for a set of parameters is typically reached 
after 20 — 50 iterations. 
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